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The atomic structure, mechanical and thermodynamic stability of vacancy clusters in Cu are 
studied by atomistic simulations. The most stable atomic configuration of small vacancy clusters 
is determined. The mechanical stability of the vacancy clusters is examined by applying uniaxial 
and volumetric tensile strain to the system. The yield stress and yield strain of the system are 
significantly reduced comparing to the prefect lattice. The dependence of vacancy formation and 
binding energy as a function of strain is explored and can be understood from the liquid-drop model. 
We find that the formation energy of the vacancy clusters decreases monotonically as a function of 
the uniaxial strain, while the formation energy increases first then decreases under the volumetric 
tensile strain. The thermodynamic stability of the vacancy clusters is analyzed by calculating the 
Gibbs free binding energy and the total probability of dissociation of the vacancy clusters at 300 
K and 900 K under uniaxial and volumetric strains. We find that although most of the vacancy 
clusters appear to be thermodynamically stable, some of the immediate sized clusters have high 
probability of dissociation into smaller clusters. 

PACS numbers: 61.72.Bb,61.72.jd,61.72. J-,61.72.-y 



I. INTRODUCTION 

Vacancies are predominant point defects in metals and 
they could have profound influence on materials proper- 
ties, ranging from mechanical strength^ 2 -* 3 -, deformation 
behavior 4,5 , kinetic transport and diffusion 6,7 , to elec- 
trical conductivity and heat capacity^, to name a few. 
Thus the study of vacancy effects on properties of tech- 
nologically important materials, such as Cu, is of great 
interest. For all these properties, stress plays crucial 
roles. For example, it is generally believed that hydro- 
static tensile stress is the driving force for the stress- void 
nucleationi^ and stress gradients are responsible for the 
void growthii. The present work touches upon the stress 
induced voiding, which is one of the most important fac- 
tors that cause failures of ultra-large integrated circuits 
with Cu-based interconnects 12 . In this case, the stress is 
built up as a result of the mismatch in thermal expansion 
coefficients between the metal lines and the surround- 
ing dielectrics^ 3 -. There are other interesting problems 
related to the void growth and coalescence under vari- 
ous loading condition s 14 i 15 i 16 i 17 i 18 i 19 . For example, ex- 
perimentally TEM observations have been performed for 
neutron- irradiated Cu at 300° C which verify the transi- 
tion from stacking fault tetrahedra (SFT) to voids, de- 
pending on the size of the vacancy cluster. These ob- 
servations are consistent to the embedded-atom-method 
(EAM) simulation s 20 ! 21 . A direct transition from a 20- 
vacancy cluster to an SFT was also found by using par- 
allel replica EAM molecular dynamics simulations 16 . We 
have recently studied void formation by using multiscale 
method and discovered a connection between vacancy- 
induced bonding cage and the tendency of void formation 
in fee metals^ 2 -. 

In this paper, we will use atomistic simulations based 
on EAM potential to (1) determine the structure of the 
most stable n- vacancy clusters where n < 21; (2) ex- 



amine the mechanical stability of Cu containing vacancy 
clusters; (3) study the dependence of vacancy-cluster for- 
mation energy and binding energy as a function of uni- 
axial and volumetric tensile stress; (4) examine the ther- 
mal stability of the vacancy-clusters at different temper- 
atures. The computational method is given in Sec II. 
The convergence of vacancy formation energy is studied 
in Sec III(A); the structure of vacancy clusters is deter- 
mined in Sec III(B) and the mechanical stability of the 
vacancy clusters is examined in Sec. III(C) for both uni- 
axial and volumetric tensile strains. Finally the thermal 
stability of the vacancy clusters is considered in III(D) 
and the conclusion is in Sec IV. 



II. COMPUTATIONAL METHOD 

A cubic supercell with periodic boundary conditions 
along three directions is used in the simulations. The su- 
percell is oriented in cubic directions, i.e., x in [100], y in 
[010] and z in [001], respectively. The atomic interaction 
is described by Ercolessi- Adams EAM potential 23 , which 
has been widely used for the metallic systems, such as 
Cuidl^d£.dld2.£2.£l£l. The atomic relaxation is carried 
out with molecular statics implemented in the LAMMPS 
codes 2 ^. It is well established that the EAM potential is 
capable to provide reasonable atomic structure and en- 
ergetics for vacancies in C u 14 i 16 i 17 . 

The vacancy cluster formation energy, which describes 
the energy cost for forming a vacancy cluster with respect 
to the perfect lattice, is computed for various vacancy 
clusters, in the presence or absence of external strains. 
Specifically, under a strain e the vacancy cluster forma- 
tion energy, E^ v e of a given n- vacancy cluster in an TV- 
atom supercell is defined as 
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where e (N — n) is the total energy of the supercell 
containing (N—n) atoms with an n- vacancy cluster under 
strain e and Eq°^(N) is the total energy of the supercell 
containing N atoms in a perfect lattice with a strain of e. 
It is convenient to use the formation energy per vacancy 
E f nv e to compare among different vacancy clusters, which 
is defined as 



all three Cartesian directions^ The dynamical matrix is 
then diagonalized and the eigen-frequencies are obtained 
to evaluate the formation entropy. The Gibbs binding 
free energy of an n- vacancy cluster at temperature T is 
determined as 
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Another important quantity calculated is the vacancy 
binding energy, E^ Y e , which describes the energy gain 
when n monovacancies are combined into an n-vacancy 
cluster. The binding energy of an n-vacancy cluster un- 
der a strain e is given by 



(3) 



A positive value of E^ v e suggests that it is energetically 
favorable for the n monovacancies to form an n-vacancy 
cluster at zero temperature under the strain. Similarly, 
one can define the binding energy per vacancy, E^ v e , to 
facilitate the comparison among different vacancy clus- 
ters, which is given by 
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In order to determine the thermodynamic stability of 
a vacancy cluster, we need to calculate dynamical matrix 
and free energy of the defect system. Under the harmonic 
approximation, the dynamical matrix under strain e is 
defined as following 
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where m is the mass of the atoms, uf and nf are the 
displacement of atom i and j in direction a and /?, re- 
spectively from its relaxed equilibrium position. E\ ot is 
the total energy of the system. The vacancy formation 
entropy is evaluated via the equatio n 26 ! 27 
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where cj^ v ' e and uo i ' e represents the harmonic vibrational 
frequency for the system with and without vacancies, 
respectively. Three acoustic phonon frequencies are re- 
moved from the sum owing to the translational invariance 
of the system. Analogous to the definition of the binding 
energy, the binding entropy is defined as 



(7) 



where Sf ve and S^ v e are the formation entropy of a 
monovacancy and an n-vacancy cluster, respectively. 

The dynamical matrix is determined by displacing the 
atoms ±0.01A from their relaxed equilibrium positions in 



The positive (negative) value of G^ v indicates the n- 
vacancy cluster is thermodynamically stable (unstable) 
against a simultaneous dissociation into n monovacan- 
cies. 



III. RESULTS AND ANALYSIS 
A. Convergence of vacancy formation energy 

A convergence test is carried out to determine the suf- 
ficient system size, and results are shown in Fig. (pQ). For 
a monovacancy, we find 6ao (ao = 3.615 A) is enough 
to reach convergency for the vacancy formation energy, 
while for 21- vacancy cluster (the largest vacancy cluster 
in this study), 20ao appears to be sufficient. Therefore, in 
all our calculations, we use 20ao as the system size. With 
a supercell of 20ao x 20ao x 20ao (32000 lattice sites), we 
found that the monovacancy formation energy is 1.2837 
eV, which agrees very well with previous calculation and 
experimen t 28 ! 29 . 
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FIG. 1: (Color online) The convergence test of the vacancy 
formation energy with respect to the system size for a mono- 
vacay (dashed line, scaled to the left axis) and a 21- vacancy 
cluster (solid line, scaled to the right axis). 



B. Atomic structure of n-vacancy clusters 

The vacancy clusters containing up to twenty-one va- 
cancies are studied in this work. For each n-vacancy clus- 



ter, several most stable configurations are considered, es- 
pecially in planar, tetrahedral, and spherical structures. 
These structures are emphasized owing to their relevance 
to vacancy loops, stacking fault tetrahedra and spherical 
voids, observed in experiment si. For each n, the struc- 
ture with the maximum binding energy is regarded as the 
most energetically stable configuration for the n- vacancy 
cluster. 

The most stable di- vacancy is the one that the two va- 
cancies are the nearest neighbors. The binding energy 
is 0.158 eV, agreeing well with the previous calculation 
and experimen t 28 ! 29 . Starting from the most stable di- 
vacancy, we can search for the most stable tri-vacancy 
by placing the third vacancy at various nearest-neighbor 
position. The configuration of the highest binding energy 
of a tri-vacancy is shown in Fig. 2(a) and E% v = 0.474 
eV. The three vacancies form an equilateral triangle ly- 
ing on {111} plane, and are mutually nearest-neighbors. 
Starting from the most stable tri-vacancy, we considered 
six most probable configurations for a four- vacancy clus- 
ter, and found that the equilateral tetrahedron formed 
by four nearest-neighbor vacancies has the highest bind- 
ing energy of Ef v =0.94 eV. The atomic structure of the 
tetra- vacancy is shown in Fig. 2(b). 

In the same manner, we explored the most stable 
atomic configurations for all vacancy clusters with n < 
21. For example, we found the most stable 5- vacancy 
cluster (penta- vacancy) is a squared pyramid with a bind- 
ing energy of Ef v = 1.372 eV, shown in Fig. 2(c). The 
most stable 6-vacancy cluster (hexa-vacancy) is an oc- 
tahedron shown in Fig. 2(d) with a binding energy of 
Eq v = 2.113 eV. Interestingly, we found that the most 
stable configuration of a 10- vacancy cluster is not an equi- 
lateral tetrahedron shown in Fig. 2(e), but rather a less 
symmetric structure shown in Fig. 2(f). For n = 13, the 
most stable structure resembles a sphere, with one va- 
cancy at the center, and 12 nearest neighbors surround- 
ing it. The radius of the sphere is y/2/2ao shown in Fig. 
2(g). For n > 13, the vacancy clusters tend to form 
3-dimensional voids, appearing more like spheres than 
tetrahedra. For example, the structure of a 19- vacancy 
cluster is shown in Fig. 2(h); it has one vacancy at center, 
and 12 nearest neighbors and 6 second- nearest neighbors 
of the central vacancy filling the sphere of a radius of a®. 

The formation energy per vacancy E^ y for all vacancy 
clusters is presented in Fig. 3. Overall, the formation en- 
ergy decreases with respect to the cluster size, which can 
be understood from the liquid-drop model for meta b 3Q i 31 . 
In this model, the void formation energy can be approx- 
imated as 
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where R ws = (3V/47T) 1 / 3 is the Wi gner-Seitz radius of 
the void with volume V : a is the surface energy density, 
and 7 is the curvature energy 30 . The volume of an n- 
vacancy cluster can be approximated as V n ~ m>o, where 
is the volume of a monovacancy. Therefor the Eq. © 



FIG. 2: The most stable atomic configurations for 3, 4, 5, 6, 
10, 13, and 19-vacancy clusters. The black circle represents 
for the vacancy and the brown circle for the atoms. Note that 
the most stable configuration for the 10- vacancy cluster is not 
the symmetrical one shown in (e) , but rather the one with less 
symmetry shown in (f). 



can be rewritten as 

Kv = °i n ~^ - C2n"*, (10) 

where c x = 47rcr(3v /47r) 2 / 3 , c 2 = 2^(3 V 47r ) 1/3 • Tak- 
ing a = 0.0932 eV/A 2 , 7 = 0.119 eV/A, and v = ajj/4, 
the Eq. (10) is plotted in Fig. ([3|), which matches very 
well the atomistic results. The value of a and 7 is taken 
from reference 3 ^, a is close to the EAM surface energy 
density of Cu: 0.08 eV/A 2 for (110) surface and 0.114 
eV/A 2 for (111) surface. Interestingly, we find that the 
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formation energy curve becomes flatter as n becomes 
larger, and approaches an asymptotic value of 0.7 eV. 
This value is very close to the monovacancy formation 
energy on a flat (111) surface of Cu, which is 0.73 eV. 
The reason for this asymptotic behavior is that as the 
vacancy cluster becomes larger, the curvature of the void 
surface becomes smaller, and thus the formation energy 
approaches the value of a monovacancy on a flat surface. 

We also present the binding energy per vacancy as a 
function of n in Fig. 3, shown in diamond curve and 
scaled to the right axis. The binding energy increases 
monotonically with n, and is positive for all vacancy 
clusters, suggesting that the monovacancies prefer to ag- 
gregate in Cu, in order to reduce the number of broken 
bonds. From Eq. ([TOj) and Eq. (|4]), we can express E^ v 
as 
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which is plotted as the dashed line in Fig. ([3]). It appears 
that the liquid-drop model reproduces the atomistic re- 
sults surprisingly well. 




FIG. 3: The formation energy (circle, scaled to left) and bind- 
ing energy (diamond, scaled to right) per vacancy of the n- 
vacancy clusters. The dot-dashed line and dashed line are 
the prediction of E f nv and from the liquid-drop model 
respectively. 



C. Mechanical stability 

In this section, we examine the mechanical stability of 
the vacancy clusters by applying external uniaxial tensile 
strain and volumetric tensile strain to the system. The 
strain-free structures of the vacancy clusters are taken 
from the results obtained in Sec. B. 



1. Uniaxial tensile strain 

The uniaxial tensile strain is applied along the y axis, 
and the x and z axes are free to change by conserving 



the volume. The strain is applied quasi- statically, with 
each loading that is 1.001 times that of the previous step. 
We examine the mechanical response of the vacancy clus- 
ters by showing the strain energy and tensile stress a yy 
as a function of applied strain e yy . Because the results 
of all vacancy clusters follow the same pattern, we select 
the 13-vacancy cluster as a representative to illustrate 
the salient features of the results. As shown in Fig. 4, 
we find that both strain energy and the tensile stress in- 
crease monotonically with strain initially, followed by a 
sudden drop at the yield point. The system behaves in a 
linear elastic fashion, i.e., the energy increases quadrat- 
ically and the stress increases linearly. There are two 
importance results: (1) the perfect and the defect lattice 
have almost identical elastic behavior - the two curves co- 
incide with each other for both energy and stress. (2) The 
presence of the vacancy defect significantly reduces the 
yield strength of the material. The yield stress and strain 
are 7.32 GPa and 0.117 for the prefect lattice, and 4.98 
GPa and 0.085 for the defect system, respectively. At 
the yield point (the drop in the curves), dislocations are 
nucleated and propagate through the system, evidenced 
by the Central Symmetry Parameter (CSP)^ 2 . plot in 

Fig. 5. CSP is defined as CSP = J2 |^ + J R i+6 | 2 /24L> 2 , 

2=1 

where Ri and Ri+Q are the vectors or bonds correspond- 
ing to the six pairs of opposite nearest neighbors in the 
deformed fee lattice 32 , and D is the nearest neighbor dis- 
tance in the perfect fee lattice. For each atom, there 
is a CSP value associated with it. Dislocation core and 
stacking fault correspond to the CSP values ranging be- 
tween 0.003 and 0.1. Just before the yielding, there is no 
dislocation in the system (Fig. 5a) while just after the 
yielding, dislocations appear (Fig. 5b). 
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FIG. 4: The strain energy (scale to the left axis) and stress 
(scale to the right axis) vs. uniaxial tensile strain for a 
13-vacancy cluster (diamond and solid line) and perfect lat- 
tice (circle and square) respectively. 

In Fig. 6, we present the yield strain (left axis) and 
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yield stress (right axis) as a function of n. Overall, the 
yield strain traces approximately the yield stress, and for 
large vacancy clusters, the yield stress/strain is signifi- 
cantly reduced comparing to the perfect lattice. How- 
ever, the vacancy effect on the mechanical strength is 
more complex than the simple monotonic behavior. More 
atomistic studies are required to unravel the complicated 
vacancy-dislocation interactions. 
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FIG. 5: CSP plot for the system immediately before (a) and 
after (b) the yield point. The color bar represents the magni- 
tude of CSP values in the system. 
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bonds are weakened owning to the stretching, therefore 
it would cost less energy to break them comparing to the 
strain- free case. The present result is consistent with that 
in Al obtained by orbital-free density functional theory 
calculations 33 . It appears that some of the energy curves 
tend to converge (or cross) at larger strains, for example, 
n = 5, 6, 10. The apparent convergence suggests that a 
large stress can change the relative stability of vacancy 
clusters of different sizes, which may trigger structural 
changes among them. 
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FIG. 7: The formation energy per vacancy under uni- 

axial stain e for several representative vacancy clusters. 

The binding energy per vacancy E^ v e is shown in 
Fig. (|HJ). The binding energy remains positive for all 
strains less than 0.062. For most of the vacancy clus- 
ters, the binding energy decreases as a function of strain 
because the attractive vacancies are pulled apart by the 
stress. It appears that large strains could lead to break- 
up or re- arrangement of the vacancy clusters. 



FIG. 6: The yield strain (circle, scaled to the left) and yield 
stress (diamond, scaled to the right) of n- vacancy clusters 
under uniaxial strain. 

We have calculated the formation energy v e of the 
vacancy clusters under uniaxial strain up to 0.062 be- 
yond which some of the defect systems start yielding. 
The results are summarized in Fig. ([7j) for representa- 
tive vacancy clusters - the results of all vacancy clusters 
have the same general features. It is found that the for- 
mation energy decreases monotonically as a function of 
strain for all clusters. This is because vacancy forma- 
tion energy results from the energy cost of bond break- 
ing around the vacancies. By applying tensile stress, the 



2. Isotropic volumetric tensile strain 

Next we study the behavior of the vacancy clusters un- 
der isotropic volumetric tensile deformations. In this case 
of triaxial loading, the volume of the system is no longer 
conserved, but with an equal strain applied along x,y,z 
directions simultaneously. At each deformation step, the 
size of the simulation box is 1.0005 times that of the pre- 
vious step. 

The strain energy and stress as a function of strain are 
shown in Fig. 9, with n = 13 as a representative because 
all vacancy clusters behave in the same way. The energy 
increases with the strain before the drop at a volume 
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FIG. 8: The binding energy per vacancy i^v,e under uniaxial 
strain e for several representative vacancy clusters 



where Vo and Kq are the volume and the bulk modulus of 
the solid at equilibrium; x = (V/V ) 1/s , £ = \{K' - 1), 



and K'q = f^|p=o- 
Vb(l + e) and .x — 1 ^ 
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For the volumetric loading, V = 
~e, thus we can rewrite Eq. ([12]) as 



K e 



(l + ^/3) : 



■exp[ 



3 J 



(13) 



Eq. (p~3|) is plotted as the triangle curve in Fig. ([9]) and 
it agrees very well with the EAM data by taking the 
experimental value of K$ — 140 GPa and £ = 5.2. 

The yield strain and yield stress for all vacancy clus- 
ters under the volumetric deformation are plotted in 
Fig. ([TO]) . The energy/stress curves of the defect systems 
trace closely to that of the perfect lattice, indicating the 
vacancy clusters do not change the behavior of the ma- 
terial before yielding. Overall, both the yield stress and 
strain decrease as a function of n, suggesting that va- 
cancy clusters can drastically reduce the yield strength 
of the material. It should be noted that the quantitative 
values of the yield stress/strain obtained here do not cor- 
respond to the realistic situation due to the high vacancy 
concentrations in the simulations; however, the qualita- 
tive features remain to be valid. 
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FIG. 9: The strain energy (scale to the left axis) and stress 
(scale to the right axis) vs. isotropic volumetric strain for a 
13- vacancy cluster (diamond and solid line) and perfect lat- 
tice(circle and square) respectively. The triangle curve repre- 
sents the EOS from Eq. (13). 



strain of e = 0.13. The drop of the energy is an indica- 
tion of the onset of plasticity of the system. The stress 
drops at the same strain as the energy. For the isotropic 



volumetric tensile strain, a y 



-P, and 



P is the pressure. The non-linear behavior of the stress- 
strain relation in Fig. 9 can be understood from the 
equation of state (EOS) of solids 34 : 

P{V) = 3Ko{ \- x) expftl-x)], (12) 
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FIG. 10: The yield strain (circle, scaled to the left) and stress 
(cross, scaled to the right) of n-vacancy clusters under the 
volumetric tensile strain. 

The formation energy per vacancy under the volumet- 
ric tensile strain is shown in Fig. ([TT]). Overall, we find 
that the formation energy first increases with the strain 
to a maximum and then decreases; this is different from 
the uniaxial case. The difference is due to the different 
loading conditions: in the uniaxial case, the volume of 
the system is fixed, therefore to a good approximation, 
the surface area of the void is also fixed. On the other 
hand, the surface area of the void increases as the vol- 
ume increases in the triaxial case, therefore the inner sur- 
face energy increases as a function of the tensile strain. 
The increasing surface energy dominates the decreasing 
bonding energy contribution at small strains, hence the 
formation energy reaches a maximum at an immediate 
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strain. For small vacancy clusters, such as n = 1 and 
n = 2, the surface energy contribution is negligible ow- 
ing to the small size of the vacancy clusters, thus their 
formation energy increases monotonically. 
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FIG. 12: The binding energy per vacancy (eV) e under 
volumetric strain e for several representative vacancy clusters. 



FIG. 11: The formation energy per vacancy i^ V)£ under the 
volumetric strain e for several representative vacancy clusters. 

For the volumetric strain, the binding energy per va- 
cancy increases monotonically with the strain, and it 
increases faster when n is larger. Such a behavior can 
be understood again from the liquid-drop model. Under 
the volumetric strain e, the volume of a monovacancy 
is vo(l + en) 3 where the linear strain en ~ e/3. Thus 
the two coefficients in Eq. (11) become ci(l + en) 2 and 
C2(l + en) respectively. E% v can then be expressed as 
a quadratic function of en or e: 
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ci(l — n 
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Since ci(l — n~z) > and [2ci(l — n~s) — 02(1 — n~s)] > 
0, ^ v e would increase monotonically with the volumet- 
ric strain e. The larger the n value, the larger the co- 
efficient of the quadratic term c\(l — n _ s) 5 hence E^ w e 
increases faster with respect to the strain. 



D. Thermodynamic Stability 

Finally, we address the thermodynamic stability of the 
vacancy clusters. At finite temperatures T, a vacancy 



cluster could dissociate into smaller clusters, or eventu- 
ally into monovacancies. In the following, we use the con- 
cept of total probability of dissociation Ptot(^) 22 to ana- 
lyze the thermal stability of the vacancy clusters. P tot (n) 
can be expressed as P tot = Pi (n) where Pi (n) repre- 
sents the probability of a particular (the zth) way of dis- 
sociation. For example, an n- vacancy cluster may disso- 
ciate into n monovacancies simultaneously, whose prob- 
ability, P\{n) is given by 



Pi(n) = 



(Civ)" 



-(nG* Vie -Gl Vte )/k B T 



J!_ e -GZ v Jk B T. 
9nv 



(15) 



or an n- vacancy cluster could break into a monovacancy 
plus an (n — l)-vacancy cluster, whose probability P 2 (^) 
is expressed as 



(14) P*[n) 



ClvC( n _l)v _ ™<7(n-l)v 



(n - l)g n 



-G 



F»-l)v,e)ABT. 



(.16) 



where g nv is the coordination number of the n-vacancy 
cluster 35 , etc. The concentration of the n-vacancy clus- 
ter is denoted by c nv and &b is the Boltzmann con- 
stant. The definition of Pi(n) can be used in both ther- 
mal equilibrium condition and a constant supersatura- 
tion condition 22 . Here we have to consider all possible 
dissociation paths in calculating a given P tot (n). The 
variation of G^ v e with respect to n determines the disso- 
ciation probability 2 ^; while an n-vacancy cluster is ther- 
mally stable under a positive and rapidly increasing G^ v e 
function, a negative and slowly increasing G^ v e would re- 
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suit in a high dissociation probability of the n-vacancy 
cluster. 

With Eqs.(5)-(7), we can evaluate the binding entropy 
Sl^v e , and then the Gibbs binding free energy of an 
n-vacancy cluster G^ v e at a given temperature T. In 
Fig. (H3D, we present G^ v e at T = 300K and T = 900K 
under the uniaxial and volumetric tensile strains. From 
the Gibbs binding free energy, one can estimate the total 
dissociation probability P to t(n) of an n-vacancy cluster. 
The calculation reveals that the total dissociation proba- 
bilities of all n-vacancy clusters, except n = 14, are much 
less than 0.5 under both strain- free and volumetric tensile 
strain conditions, which means they are thermodynam- 
ically stable. On the other hand, under the volumetric 
strain, the 14- vacancy cluster has more than 50% proba- 
bility to dissociate into a 13- vacancy cluster and a mono- 
vacancy because G? 4v is smaller than Gf 3v e . Under 
a uniaxial tensile strain, 10- and 11-vacancy clusters are 
not stable at T = 900K; the 10- or 11-vacancy cluster has 
a high probability to dissociate into two 5-vacancy clus- 
ters or one 5- and one 6- vacancy clusters, respectively. 
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FIG. 13: The Gibbs binding free energy of n-vacancy clusters 
at different temperatures for strain- free (circle), uniaxial ten- 
sile strain (diamond) and volumetric tensile strain (square). 



IV. CONCLUSION 

We have examined structure, mechanical and thermo- 
dynamic stability of n-vacancy clusters where n goes from 
2 to 21. We determine the most stable atomic structure 
of the vacancy clusters based on simple intuitive con- 
siderations and energetic calculations. More extensive 
(such as global search) and accurate (such as quantum 
mechanical) calculations are required to settle the issue 
with more confidence. The present work should only be 
considered as a preliminary effort. We find that the for- 
mation energy per vacancy decreases while the binding 
energy increases as a function of cluster size n. We probe 
the mechanical stability of the vacancy clusters by apply- 
ing uniaxial and volumetric tensile strains. It is observed 
that yield stress and yield strain of the material are sig- 
nificantly reduced by the vacancy clusters. However, the 
presence of the vacancy cluster does not change the elas- 
tic behavior of the material before yielding. We find that 
formation energy per vacancy decreases as a function of 
uniaxial strain for all clusters, while the binding energy 
shows a more complicated behavior. For volumetric de- 
formations, the formation energy per vacancy increases 
first then decreases as a function of strain. The increase 
of the formation energy at smaller strains is due to the 
increased surface energy associated with the void. We 
determine the thermodynamic stability of the clusters by 
calculating the Gibbs free binding energy and resultant 
probability of dissociation. We find that most of the va- 
cancy clusters under study are thermodynamically stable 
except the 14- vacancy cluster, which has a high probabil- 
ity of dissociating into a 13- vacancy cluster and a mono- 
vacancy. In addition, at 900K the 10-and 11-vacancy 
clusters will dissociate into two 5-vacancy clusters or a 5- 
and a 6-vacancy clusters, respectively under a uniaxial 
strain. 
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